Final Fate of Subcritical Evolutions of Boson Stars 
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We present results from a study of Type I critical phenomena in the dynamics of general rela- 
tivistic boson stars in spherical symmetry. The boson stars are modelled with a minimally coupled, 
massive complex field (with no explicit self-interaction), and are driven to the threshold of black 
hole formation via their gravitational interaction with an initially implodingpulse of massless scalar 
field. Using a distinct coordinate system, we reproduce previous results [2,12] 5 including the scaling 
of the lifetime of near-critical configurations, as well as the fact that such configurations are well 
described as perturbed, one-mode-unstable boson stars. In addition, we make a detailed study of 
the long-time evolution of marginally subcritical configurations. Contrary to previous claims [f|, Q], 
we find that the end state in such cases does not involve dispersal of the bulk of the boson star 
field to large radial distances, but instead can be generically described by a stable boson star ex- 
ecuting large amplitude oscillations. Furthermore we show that these oscillations can be largely 
identified as excitations of the fundamental mode associated with the final boson star, as computed 
in perturbation theory. 
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I. INTRODUCTION 

Over the past decade or so, intricate and unexpected 
phenomena related to black holes have been discovered 
through the detailed numerical study of various mod- 
els for gravitational collapse, starting with one of the 
authors' investigation of the spherically symmetric col- 
lapse of a massless scalar field These studies gen- 
erally concern the threshold of black hole formation (a 
concept described below), and the phenomena observed 
near threshold are collectively called (black hole) critical 
phenomena, since they share many of the features asso- 
ciated with critical phenomena in statistical mechanical 
systems. The study of critical phenomena continues to be 
an active area of research in numerical relativity, and we 
refer the interested reader to the review article by Gund- 
lach [1] for full details on the subject. Here we will simply 
summarize some key points that are most germane to the 
work described in this paper. 

To understand black hole critical phenomena, one must 
understand the notion of the "threshold of black hole for- 
mation" . The basic idea is to consider families of solu- 
tions of the coupled dynamical equations for the grav- 
itational field and the matter field that is undergoing 
collapse (a complex scalar field, <p, in our case). Since 
we are considering a dynamical problem, and since we 
assume that the overall dynamics is uniquely determined 
by the initial conditions, we can view the families as be- 
ing parametrized by the initial conditions — variations in 
one or more of the parameters that fix the initial values 
will then generate various solution families. We also em- 
phasize that we are considering collapse problems. This 
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means that we will generically be studying the dynamics 
of systems that have length scales comparable to their 
Schwarzschild radii, for at least some period of time dur- 
ing the dynamical evolution. We also note that we will 
often take advantage of the complete freedom we have 
as numerical experimentalists to choose initial conditions 
that lead to collapse, but which may be highly unlikely 
to occur in an astrophysical setting. 

We now focus attention on single parameter families of 
data, so that the specification of the initial data is fixed 
up to the value of the family parameter, p. We will gener- 
ally view p as a non-linear control parameter that will be 
used to govern how strong the gravitational field becomes 
in the subsequent evolution of the initial data, and in par- 
ticular, whether a black hole forms or not. Specifically, 
we will always demand that any one-parameter family of 
solutions has the following properties: 

1. For sufficiently small values of p the dynamics re- 
main regular for all time, and no black hole forms. 

2. For sufficiently large values of p, complete gravi- 
tational collapse sets in at some point during the 
dynamical development of the initial data, and a 
black hole forms. 

From the point of view of simulation, it turns out to be 
a relatively easy task for many models of collapse to con- 
struct such families, and then to identify two specific pa- 
rameter values, p~ (p + ) which do not (do) lead to black 
hole formation. Once such a "bracket" [p~,p + ] has been 
found, it is straightforward in principle to use a technique 
such as binary search to hone in on a critical parameter 
value, p* , such that all solutions with p < p* (p > p*) do 
not (do) contain black holes. A solution corresponding 
to p = p* thus sits at the threshold of black hole forma- 
tion, and is known as a critical solution. It should be 
emphasized that underlying the existence of critical so- 
lutions are the facts that (1) the end states (infinite-time 
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behaviour) corresponding to properties 1 and 2 above are 
distinct (a spacetime containing a black hole vs a space- 
time not containing a black hole) and (2) the process 
characterizing the black hole threshold (i.e. gravitational 
collapse) is unstable. We also note that we will term evo- 
lutions with p < p* subcritical, while those with p > p* 
will be called supercritical. 

Having discussed the basic concepts underlying black 
hole critical phenomena, we now briefly describe the fea- 
tures of critical collapse that are most relevant to the 
research described below. 

First, critical solutions do exist for all matter models 
that have been studied to date, and for any given matter 
model, almost certainly constitute discrete sets. In fact, 
for some models, there may be only one critical solution, 
and we therefore have a form of universality. 

Second, critical solutions tend to have additional sym- 
metry beyond that which has been adopted in the speci- 
fication of the model (e.g. we will impose spherical sym- 
metry in our calculations). 

Third, the critical solutions known thus far, and the 
black hole thresholds associated with them, come in two 
broad classes. The first, dubbed Type I, is characterized 
by static or periodic critical solutions (i.e. the additional 
symmetry is a continuous or discrete time-translational 
symmetry), and by the fact that the black hole mass just 
above threshold is finite (i.e. so that there is a minimum 
black hole mass that can be formed from the collapse). 
The second class, called Type II, is characterized by con- 
tinuously or discretely self-similar critical solutions (i.e. 
the additional symmetry is a continuous or discrete scal- 
ing symmetry) , and by the fact that the black hole mass 
just above threshold is infinitesimal (i.e. so that there is 
no minimum for the black hole mass that can be formed). 
The nomenclature Type I and Type II is by analogy with 
first and second order phase transitions in statistical me- 
chanics, with the black hole mass viewed as an order 
parameter. 

Fourth, solutions close to criticality exhibit various 
scaling laws. For example, in the case of Type I col- 
lapse, where the critical solution is an unstable, time- 
independent (or periodic) compact object, the amount 
of time, r, that the dynamically evolved configuration is 
well approximated by the critical solution per se satisfies 
a scaling law of the form 

r(p) 7ln|p-p*| , (1) 

where 7 is a universal exponent in the sense of not de- 
pending on which particular family of initial data is used 
to generate the critical solution, and ~ indicates that the 
relation (JTJ) is expected to hold in the limit p — ► p* . 

Fifth, and finally, much insight into critical phenomena 
comes from the observation that although unstable, crit- 
ical solutions tend to be minimally unstable, in the sense 
that they tend to have only a few, and perhaps only one, 
unstable modes in perturbation theory. In fact, if one 
assumes that a Type I solution, for example, has only a 



single unstable mode, then the growth factor (Lyapunov 
exponent) associated with that mode can be immediately 
related to the scaling exponent 7 defined by (Q]). 

In this paper we will be exclusively concerned with 
Type I critical phenomena, where the threshold solutions 
will generally turn out to be unstable boson stars. Pre- 
vious work relevant to ours includes studies by Hawley 
Q and Hawley & Choptuik of boson stars in spher- 
ically symmetry. We extend this work and show that, 
contrary to previous claims [l], 0] that subcritical solu- 
tions disperse most of the original mass of the boson star 
to large distances — the late time behaviour of subcritical 
evolution is characterized by oscillation about a stable 
boson star solution. We also apply a linear perturba- 
tion analysis similar to that in [2| and confirm that 
the observed oscillation modes agree with the fundamen- 
tal modes given by perturbation theory. (We use a code 
kindly provided by S. Hawley [5[ to generate the frequen- 
cies from the perturbation analysis.) 

The outline of the rest of this paper is as follows: 
111 Section Hfl we describe the mathematical formulation 
for our numerical simulations, which includes III Al the 
model for the boson stars, and lll Bt the initial value prob- 
lem. In Section [IIII we present results of our simulations: 
in IIII Al we present the setup of numerical experiments, 
in IIII Bl the Type I character of the critical solutions is 
demonstrated. [HI Cl contains a discussion of the end state 
of subcritical evolutions and is followed by some pertur- 
bation analysis in llll Dl Section ITVl summarizes our find- 
ings, while the finite difference approximations used and 
our convergence testing of our implementations of them, 
are given in Appendices O and [B] respectively. 

In what follows we base our work in the context of clas- 
sical field theory, and we choose units in which G = c = 1 . 
In addition, without loss of generality, we restrict our- 
selves to the case where the particle mass, m, associated 
with the complex scalar field satisfies m = 1. 



II. MATHEMATICAL FORMULATION 

A. The model 

Our model for boson stars involves a self-gravitating 
massive complex scalar field, <j) = 4>i + i(f>2, minimally 
coupled to gravity as given by general relativity. (Note 
that we do not make the complex scalar field explicitly 
self-interacting — in the literature, the stationary config- 
urations in this case are sometimes called "mini" bo- 
son stars.) An additional, massless real scalar field, </>3, 
also minimally coupled to gravity, is used to dynamically 
"perturb" the boson star. The interaction between the 
massive complex scalar field and the massless real scalar 
field is thus through the gravitational field alone. The 
whole system can be described by the action 
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(WV M 0* + m 2 (jxt)*) 
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where R is the spacetime Ricci scalar and m is the mass 
of the bosonic particle. Variations of the action with re- 
spect to the metric, g M „, the complex scalar field, 4>, and 
the real scalar field, 03, yield the equations of motion, 
which are the Einstein equation, the Klein-Gordon equa- 
tion and the wave equation, respectively: 



ds 2 = {-a 2 + i/> 4 /3 2 ) dt 2 + 2i/)*/3dtdr 

+ t// (dr 2 + r 2 dVL 2 ) , (9) 

where a, j3 and ip are the lapse function, r-component 
of the shift vector and the conformal factor respectively, 
and all are functions of t and r. We further define new 
variables to transform the Klein Gordon and wave equa- 
tions into a first order (in time) system: 



(10) 



Rfj,u — -g^vR — SttT^l, , 



(3) 



n, 



->i' 2 



(4>i - /30-) , 



(11) 



and 



where 



VV^ - m 2 (t) = , 



VV^ = 0. 



rp _ rp(j> i rp(p 3 



(4) 



(5) 



(6) 



where i = 1, 2 or 3, ' = d/dr and ' = d/dt. As with the 
geometric variables, 4>i. $i and IT are functions oft and 
r alone. 

With these definitions, the Hamiltonian constraint and 
momentum constraints are given by Q 



A A ( r ^\ + ±K r r 2 = 
i/; 5 dr 3 \ dr J 16 r 



' Hi (^ + n?) 

^ 4 



(12) 



-<^( V Q 0V Q r +™ 2 |0| 2 )] , (7) 



and the Klein-Gordon and wave equations become 



(13) 



-9nv V a </> 3 V Q <; 



(8) 



Equations ©-([HI) completely determine the dynamics 
of our system (up to coordinate transformations), once 
appropriate initial conditions and boundary conditions 
are specified. 

To study the system numerically, we adopt the stan- 
dard "3+1" ADM formalism [1,0.* Since we restrict our- 
selves to spherically symmetry the metric can be written 
in a much simpler form than in the generic case. Here 
we use maximal-isotropic coordinates, which is a differ- 
ent system than that used in P, 0] ■ We note in passing 
that although the accuracy of finite difference calcula- 
tions in any given coordinate system can in principle be 
estimated using intrinsic means (e.g. convergence tests), 
we feel that it is nonetheless useful to reproduce the cal- 
culations of [l|, in a distinct coordinate system. 

In maximal-isotropic coordinates, the line element can 
be written as: 



:IT + , 



n, 



3 d 



aK r „ 



r^ 4 1/311* + — $ 



(14) 
(15) 

aip 2 m 2 4> % (1 - <5 l3 ) 
(16) 



In addition to equations (|12|) - p6)) . we need to de- 
termine the lapse function and shift component using 
our specific coordinate choices. The maximal condition, 
which maximizes the 3- volume of the t = const, slices, is 
given by K = K l j = 0. This is implemented by choosing 
initial data so that K(0, r) = and then demanding that 



K(t,r) = 0. 



(17) 



for all t and r. 
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This leads to the following linear ODE for a(t,r), 
which must be solved at each instant in time in order 
to maintain the maximal condition [8J|: 



"" + -^-7^7 (r 2 ip 2 ) a' 



rip 2 dr 2 



4™v 4 E # - 8 - E n ' - |(^ r « 



0. 



(18) 



The isotropic condition, which is implicit in the cho- 
sen form of the metric ([9]) , demands that the 3- metric of 
each t = const, hypersurface be conformally flat. This re- 
quirement leads to the following ODE for the shift vector 
component, (3{t, r): 



-aK r , 



(19) 



Equations Q12 p -(|19 p constitute a complete set of dif- 
ferential equations governing our model. Note that our 
approach is an instance of so-called fully constrained evo- 
lution, wherein all of the geometric variables — ip and K r r 
in this case — are computed at each time step using con- 
straint equations. To completely fix a solution — given 
initial data — we also must impose regularity and bound- 
ary conditions at r = and r — ► oo respectively. 

Regularity at the origin, r = 0, requires 



V'(t,0) 
K r r (t,0) 
<*'(*, 0) 

#(*>o) 
nj(t,o) 



o. 

o, 
o, 
o, 

0. 



(20) 
(21) 
(22) 
(23) 
(24) 



4\ 



n, 



o. 



n< + nj + — = o, 



(28) 
(29) 



for some functions C(t) and -D(i). The last two of these 
equations are approximate Sommcrfcld conditions that 
assume that as r — > oo, the three scalar field compo- 
nents, <pi, are purely outgoing with amplitudes decaying 
as 1/r. For given initial data, eqs. (fT2"|) - (f!?9")) now com- 
pletely determine our system. 

For diagnostic purposes, we also define the mass aspect 
function 



M(t, r) 



ip 2 r 



K r „ 



2iP'r 2 (ip + nf)') , (30) 



which is equal to the ADM mass in any vacuum region 
exterior to the support of matter. 

In addition, although ip and K r r are ultimately deter- 
mined from the constraint equations, the following evo- 
lution equations are used for providing initial estimates 
for the iterative constraint-solving process. 



1 

V 2 v r 2ip 



(31) 



K\ = 0K r ' r 



2a 



(rip 2 )' 2 r 2 ip e 
8nm 2 a\(j)\ 2 . 



or 



(32) 



Full details of our finite differencing scheme are given 
in App. El 



whereas the outer boundary conditions are 



B. The initial value problem 



lim ip(t, r) = 1 

r — >oo 

lim a(t, r) = lim 

r — >oo r — >oc 

= 1 - 



C(t) 



+ 0(r- 2 ) 



> ip(t, r 
2C(t) 



(25) 



+ 0(r~ 2 ), (26) 



lim (3{t,r) = ^ + 0(r- 2 ), (27) 



and 



The primitive object in our model problem is a 
(ground-state) boson star, represented by a configura- 
tion of massive complex scalar field, centred at the origin. 
Ideally one would like a "star" to be described by a lo- 
calized, time-independent matter source that generates 
an everywhere regular (i.e. non-singular) gravitational 
field. However, for the case of a complex scalar field, it 
can be shown that such regular, time-independent con- 
figurations do not exist [9J. Despite this fact, since the 
stress-energy tensor depends only on the modulus of 
the scalar field (and the gradients of the modulus), one 
can construct scalar field configurations with harmonic 
time-dependence that produce time-independent metrics. 
Specifically, we adopt the following ansatz for boson stars 
in spherical symmetry: 
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while at the outer boundary, we have 



<f>(t,r) = &o{r)e- 



(33) 



and then demand that the spacetime be static, i.e. we 
demand that the metric admits a timclikc Killing vector 
field, x, which is orthogonal to the t = const, surfaces. 
Adapting coordinate time to the timelike Killing vector 
field, we have 



(3 = 0, 



(34) 



for all time t. Additionally, we have that the time deriva- 
tives of any of the geometrical variables identically van- 
ish. It then follows immediately that 8] 



(35) 



K \ = . 



As is necessary for the consistency of the ansatz (1331) . 
the isotropic condition for j3 (fT§]) is automatically satis- 
fied, and we are left with geometrical variables a(0,r), 
ip(0, r) and 4>o(r) that need to be determined from the 
maximal slicing condition (|18p . the Hamiltonian con- 
straint (fT2")l and the Klein-Gordon equation (|TS)) . respec- 
tively: 



2* 



■0$ [ —z+m 



2 I J.2 



= 

= A, 
= -2 



2 A 2* 
- + - + — 

r a ip 



- 4- ) <T> - V 4 ( m 2 - ^ 



(36) 
(37) 
(38) 
(39) 
(40) 



{I + |) A + 4 ^ 4 " " ^ (41) 



Here, in order to simplify notation, we have dropped 
the subscript "0" , making the identifications 0(r) = 
4>o(r) and <I>(r) = 0'(r) = <fi' (r). We have also intro- 
duced auxiliary variables \&(r) = ip'(r) and A(r) = a'(r) 
in order to cast the above system of nonlinear ODEs in a 
canonical first-order form. We assert that for any given 
value of (f)(0) = 0o(O), the system ([36]) -(|4i" j) constitutes an 
eigenvalue problem with eigenvalue lo = w(0(O)). That 
is, for any specific value of 0(0) (which one can loosely 
view as being related to the central density of the star), 
a solution of (|33|) that satisfies the appropriate regularity 
and boundary conditions will only exist for some specific 
value of lo. The system (|3l)l) - (|4"l"j) must be supplemented 
by boundary conditions, some of which are naturally ap- 
plied at r = 0, with the rest naturally set at r = oo. In 
particular, regularity at r — implies 



*(0) = 0. 
*(0) = o, 
A(0) = 0, 



(42) 
(43) 
(44) 



lim ip(r) 

r — >oo 

lim (f>(r) 

r — >oo 

lim a(r) 



C 



0. 



1-1. 



(45) 
(46) 

(47) 



Here the second condition follows from the expectation 
that 4> should decay exponentially [10( as r — > 0. 

We further note that due to the homogeneity and lin- 
earity of the slicing equation, we can always arbitrarily 
(and conveniently) choose the central value of the lapse 
via 



a(Q) 



1 



(48) 



and then, after integration of (|36[) - (l4"T]) , can rescale a 
and lo simultaneously to satisfy the true outer boundary 
condition, (I47p . for a: 



a{r) 
io(r) 



ca(r) . 
clo(t) . 



where c is given by 



2/V>(r„ 



a(r max ) 



(49) 
(50) 



(51) 



and 



is the radial coordinate of the outer boundary 



of the computational domain. 

As mentioned above, any solution of ((3"6")) - (|4"l"j) can be 
conveniently labelled by the central value of the modulus 
of the scalar field, <f>o(0) = 0(0). For any given value 
of 0o(O), we must then determine the eigenvalue, lo, and 
in the current case of maximal-isotropic coordinates, the 
central value of the conformal factor "0(0), so that all of 
the boundary conditions are satisfied. In principle, we 
can compute pairs [w,-0(O)] as a function of 0n(O) using 
a two-parameter "snooting" technique [ll], [13] ■ 

Alternatively, in some cases we generate boson star ini- 
tial data in maximal-isotropic coordinates by first con- 
structing the stars in so-called polar-areal coordinates, 
and then performing a coordinate transformation on the 
resulting solution. 

Polar-areal coordinates, which have seen widespread 
use in spherically symmetric computations in numerical 
relativity, can be viewed as the generalization of the usual 
Schwarzschild coordinates to time- dependent, spherically 
symmetric spacetimes. As with maximal slicing, the slic- 
ing condition in this case — known as polar slicing — is ex- 
pressed as a condition on the mean extrinsic curvature: 



K = K T r . 

Since in general we have K — K l i = K r , 
condition is implemented by requiring 

K 6 g(t,r) = k\(t,r)=0, 



(52) 
2K%, this 

(53) 
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for all t and r. 

The spatial coordinates are fixed by demanding that 
the coordinate r measure proper surface area (i.e. that it 
be an areal coordinate). It can be shown that this choice 
of r, together with polar slicing, further imply that (3 = 0, 
so that the line element becomes 



-a 2 dt 2 + a 2 dr 2 



'-dn 1 



(54) 



As before, to construct star-like solutions, we adopt the 
time-harmonic ansatz (|33p for the complex scalar field, 
adapt the time coordinate to the timelike Killing vector 
field, and require the spacetime to be static. We again 
find that the extrinsic curvature tensor vanishes identi- 
cally (so that, for static data, the slicing is polar as well 
as maximal), and that the momentum constraint Q13[) is 
automatically satisfied. 

Again, considering the Hamiltonian constraint, the 
Klein-Gordon equation, and the slicing condition 

K% = , (55) 
at t = 0, we have (dropping the subscript 0's as before): 



2 2 2 

a m 




(56) 



- (l + a 2 - Aitr 2 a 2 m 2 (h T 



4> 



— — m 4>a 



In this case, the regularity conditions are 

o(0) = 1, 
*(0) = 0, 

while the outer boundary conditions are 

lim 4>(r) ~ , 

r — >oo 

lim a(r) 



1 



a(r) 



(57) 
(58) 

(59) 



(60) 
(61) 



(62) 
(63) 



As before, we can convert the last condition to an inner 
condition on a by taking advantage of the linearity and 
homogeneity of the slicing equation. Specifically, we can 
again choose a(0) = 1, and then after integration of (|56p- 
(159")) simultaneously rescale a(r) as well as the eigenvalue, 
u>, so that (1631) is satisfied. 



We again consider the family of boson star solutions 
parametrized by the central value of the modulus of the 
scalar field, ^o(O). In this case, given a value of </>o(0), 
and using the conditions a(0) = 1, a(0) = 1, $(0) = 
0, we need only adjust the eigenvalue to itself in order 
to generate a solution with the appropriate asymptotic 
behaviour (i.e. so that linv— >oo <f>(r) = 0). This is a classic 
Tparameter shooting problem, which is comparatively 
easier than the 2-parameter shooting method described 
above. 

Once we have computed a solution in areal coordinates, 
we can perform a coordinate transformation from areal 
coordinates to isotropic coordinates [HI, [3] (recall that 
the maximal and polar slices coincide for the static case). 
Essentially this amounts to solving an ODE of the form 



\R=R„ 



R=R„ 



dr 
dR 



R 



(64) 



We emphasize that (1551) - (|4"Tj) or (|55|) - ([591) are used for 
generating initial data describing a spacetime that has 
no matter content other than a single boson star. To 
"perturb" a given boson star, and, in particular, to drive 
the star to the threshold of black hole formation, we im- 
plode a (spherical) shell of massless scalar field onto it. 
Specifically, we choose initial data for the massless field 
of the following "gaussian" form: 



h(0,r) = A 3 exp 



r-r 



(65) 



where A^,rQ and a are adjustable parameters, control- 
ling the overall amplitude, position and width, respec- 
tively, of the imploding gaussian wave packet. To en- 
sure that the massless field is almost purely in-going 
at the initial time, we specify the "conjugate" variable 



n. 



ip 2 /a U> 3 



3 — /3<j) 3 ' ) as follows: 



n 3 (0,r) = - $ 3 (0,r) + 



(66) 



In all of our studies described below, we have fixed tq 
and a in ([65]) to tq — 40 and a = 5. This ensures that 
the support of the massless field is well separated from 
that of the complex field (i.e. from any of the boson stars 
per se that we study) at the initial time. 

Once the complex scalar field, <j>, and the real scalar 
field, ^3, are known, the initial data for the functions 
ip(0,r), K r r (0,r), a(0,r) and /?(0,r) are computed by 
solving the Hamiltonian constraint (fT2)) . the momen- 
tum constraint (|13p . the slicing condition (fT5)) . and the 
isotropic condition respectively. 
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III. RESULTS 
A. Setup of numerical experiments 

The PDEs solved in the simulations discussed here are 
those listed in the previous section. We also provide a 
summary of the equations of motion of the system, the 
boundary conditions, and details of the finite difference 
approximation used in App. [SJ Additionally, results of 
convergence tests of our code are discussed in App. [SJ 

In order to study critical behaviour in the model we 
start with initial data for the complex field that repre- 
sents a boson star on the stable branch (i.e. a star with 
a central scalar field value 0o(O) < 0.08, using our units 
and conventions). We generally choose a configuration 
that is reasonably relativistic, i.e. with 0o(O) bounded 
away from 0, but not too close to the instability point, 
0(0) « 0.08. 

Details of the initial data setup were described in the 
previous section. Typical evolution of such initial data 
proceeds as follows. Once we have fixed the boson star 
configuration, we complete the specification of the mass- 
less scalar field initial data by fixing the overall ampli- 
tude factor, A3, and then evolve the system. Initially, 
the shell of massless scalar field implodes towards r = 
at the speed of light, while the boson star "sits" in its 
static state centered at the origin. As the in-going mass- 
less shell reaches the region of space occupied by the 
boson star, its contribution to the overall gravitational 
field tends to compress the boson star to a higher mean 
density and smaller radius. The massless field passes 
through the origin and then "explodes" outward, even- 
tually propagating off the computational domain. De- 
pending on the strength of the perturbation from the 
massless field, we find that the compressed boson star ei- 
ther relaxes to something resembling a stable boson star 
with large-amplitude oscillations, or collapses to form a 
black hole. Thus by adjusting the massless scalar am- 
plitude factor, A3 — which we generically use as the ad- 
justable parameter, p, in our study of critical behaviour 
in the model — we can tune the evolution to the thresh- 
old of black hole formation. In practice we use a bisec- 
tion search to refine our estimate of the critical value, 
A3, and can carry the search to machine precision, so 
that AA3/A3 ~ 10~ 15 using standard 8-byte floating- 
point arithmetic. Unless otherwise specified, the com- 
putations described below have been performed with a 
spatial mesh spacing Ar = 50/1024 « 0.049, a Courant 
factor At/ Ar — 0.3, and the coefficient of Kreiss-Oliger 
dissipation, td = 0.5 (see App. lAlfor the definition of 64). 

We note that our numerical calculations generate en- 
tire families of critical solutions, fundamentally reflecting 
the fact that there is a continuum of one-mode unstable 
boson star configurations (see Fig. [3]). In addition, for 
any fixed initial boson star state, the specifics of the ob- 
served threshold solution will depend on the details of 
the "perturbing" scalar field. This last fact is, however, 
irrelevant to the conclusions that we draw from our study. 



In the following section we discuss results from detailed 
studies of black hole threshold solutions generated from 
several distinct initial boson star states. Table I summa- 
rizes the values of 0o(O) that were used, the approximate 
values of A3 required to generate a critical solution, the 
location, r max , of the outer boundary of the computa- 
tional domain, and the figures that display results asso- 
ciated with the respective calculations. Since we will not 
dwell on this point below, we note that all of our calcula- 
tions confirm the basic picture previously reported that 
the black holes that form just above threshold in this 
type of collapse generically have finite mass (i.e. that the 
critical transition is Type I). 



B. Critical phenomena 

We start by examining results from a critically per- 
turbed boson star having an unperturbed central field 
value </>o (0) = 0.05. As just described, the critical mass- 
less amplitude factor, A3 ~ 0.0032, was determined by 
performing a bisection search on A3, to roughly machine 
precision. (Recall that each iteration in this search in- 
volves the solution of the time-dependent PDEs for the 
model for a specific value of A3, with all other parame- 
ters held fixed, and the criterion by which we adjust the 
bisection bracket is whether or not the simulation results 
in black hole formation.) 

A series of snapshots of dM(t,r)/dr (where M(t,r) 
is the mass aspect function) for a marginally subcritical 
evolution is shown in Fig.[TJ Full analysis of the results of 
this simulation indicate that the boson star enters what 
we identify as the critical state at t » 130, and remains in 
that state until t s» 510. It is worth noting that the boson 
star actually completes its collapse into a more compact 
configuration well after the real scalar field has dispersed 
from the boson star region. We also note that the amount 
of time, r, spent in the critical state — r ~ 380 in this 
case — is a function of how closely the control parameter 
has been tuned to criticality. Specifically, we expect r 
to be linear in In | A3 — A3 1 (see ([T])), and we will display 
evidence for this type of scaling below. 

Fig. [U shows the time evolution of the central modu- 
lus of the complex scalar field for marginally subcritical 
evolutions generated from boson star initial states with 
0o (0) = 0.035, 0.04 and 0.05. From the figure we can see 
that in all three cases the perturbed stars enter an ex- 
cited, critical state at t w 100 and remain in that state for 
a finite time which is a function of 0o(O) (i.e. of the initial 
state). Additionally, at least for the cases 0o(O) = 0.035, 
0o (0) = 0.04, the figure provides evidence that follow- 
ing the critical evolution phase, the excited stars relax to 
states characterized by large amplitude oscillations of the 
complex field. This behaviour will be examined in more 
detail below. Finally, also apparent in the plot are the 
smaller-amplitude oscillations that occur during the pe- 
riods of critical evolution. Previous work [l], 0] indicated 
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TABLE I: Summary of parameters used to generate the results displayed in Figs. IH8I Listed for each distinct computation or 
numerical experiment are the relevant figure numbers, central amplitude of the complex field, (f>o{0), the overall massless scalar 
amplitude factor, A% (see I65p . that generates a marginally-critical solution, and the maximum radial coordinate, r max , of the 
computational domain. Other parameters defining the massless scalar initial profile (|65[1 are held fixed at ro = 40, a — 5 for 
all simulations. Other numerical parameters are chosen to be Ar = 50/1024 ps 0.049, At/ Ar = 0.3 and e<j = 0.5, and are also 
fixed for the calculations discussed here. 



that these oscillations can be interpreted as excitations 
of the (stable) first harmonic mode of the unstable boson 
star that is acting as the critical solution — the unstable 
fundamental mode is the one that determines whether 
or not the configuration will evolve to a black hole. Al- 
though we have not studied this matter in any detail, 
we assume that the same picture holds for our current 
calculations. 

The results from our simulations of critically perturbed 
boson stars are thus in agreement with the previous stud- 
ies [l|,[2| which identified the critical states as excited (pri- 
marily in the first harmonic mode) , unstable boson stars. 
Following that work we display in Fig. [3] an approximate 
correspondence between the initial boson stars and the 
critical solutions. The solid line traces the one-parameter 
family of static boson stars (parameterized as usual by 
</>o(0)), where we have defined the radius, R, of a boson 
star so that M(R) = 0.99 M(oo) = 0.99 M ADM . The 
triangles indicate the initial stable boson star configura- 
tions, the squares indicate our best estimate of the cor- 
responding unstable critical boson star states, and each 
arrow schematically depicts the transition between the 
two states that is induced by the perturbing scalar field. 
We note that to identify which unstable boson star is act- 
ing as the critical solution — which is equivalent to iden- 
tifying an effective value of <j>o(0) — we time average the 
central modulus of the complex field, \<fc(t, 0)| during the 
period of critical evolution. In addition, in accord with 
previous results, we observe that in all cases the mass of 
the unstable critical state is larger than that of the pro- 
genitor boson star, indicating that a significant amount 
of mass-energy is extracted from the massless scalar field 
through its purely gravitational interaction with the com- 
plex field. 

As discussed previously, for both subcritical and su- 
percritical simulations, the closer one tunes ^3 to the 
critical value A\, the longer the perturbed star will per- 
sist in the critical state. Specifically, we observe scaling 
of the lifetime, r, of the critical evolution of the form 

t(A 3 )~-j\u\A 3 -A* 3 \, (67) 



where we define the lifetime to be the lapse of coordinate 
time from the start of the evolution, t = 0, to the time of 
first detection of an apparent horizon, and where 7 is a 
scaling exponent that depends on which of the infinitely 
many one-mode unstable boson stars acts as the criti- 
cal solution in the particular scenario being simulated. 
We note that the details of the definition of r are not 
important to the determination of 7 in (|67p since 7 ac- 
tually measures the differential in lifetime with respect 
to changes in A 3 — A 3 , and this differential is insensitive 
to precisely how we define t, at least as ^3 — * A\. In 
addition, we note that in using coordinate time in our 
definition of the scaling relationship (|67[) , we are defining 
the scaling with respect to proper time at spatial infin- 
ity. Another choice — arguably more natural — would be 
to define r in terms of the proper time measured by an 
observer at rest at r = (central proper time). Since 
the critical solutions are nearly static, the relation be- 
tween these two different definitions of time would be a 
specific factor for each distinct value of 4>o(0), and would 
thus lead to a <f>$ (O)-dependent "renormalization" of the 
scaling exponents, 7. 

Fig. 2] shows measured scaling laws from supercritical 
evolutions of perturbed boson stars defined by 4>o(0) = 
0.02,0.035,0.04 and 0.05. It is clear from these plots 
that, at least as A 3 — > A 3 , we have lifetime scaling 
of the form (|6"T)) . Estimated values of 7 — computed 
from linear least-squares fits to the plotted data — are 
7 = 8.1,11,14,17 for 0o (0) = 0.02,0.035,0.04,0.05, re- 
spectively. We note that according to the now stan- 
dard picture of critical collapse (see for example [H), 
each value of 7 can be identified with the reciprocal Lya- 
punov exponent (i.e. growth factors) of the single unsta- 
ble mode associated with the corresponding critical solu- 
tion. Again, the reason that we observe different values 
of 7 for different choices of initial boson star (different 
values of 4>o(0)) is that distinct critical solutions are be- 
ing generated in the various cases. That is, we cannot 
expect universality (with respect to initial data) in this 
case because the model admits an entire family of one- 
mode unstable solutions that sit at the threshold of black 
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FIG. 1: Critical evolution of a perturbed boson star with 
<f>o (0) = 0.05 and mass, Madm = 0.62 (using our units and 
conventions) . This figure shows the time development of con- 
tributions to dM/dr from the complex (solid line) and real 
(dashed line) scalar fields. Note that the temporal spacing 
between successive snapshots is not constant — the time in- 
stants displayed have been chosen to illustrate the key fea- 
tures of the near-critical evolution. Also note that we have 
multiplied the value of dM/dr for the real scalar field by a 
factor of 8 to aid in the visualization of that field's dynam- 
ics. The evolution begins with a stable boson star centered 
at the origin, and an in-going gaussian pulse (shell) of mass- 
less, real scalar field that is used to perturb the star. The 
overall amplitude factor, A3, of the initial real scalar field 
profile (see Q65p). is the control parameter for generating the 
one-parameter family of solutions that interpolates through 
the black hole threshold. For the calculation shown here, A3 
has been tuned to a critical value A3 » 0.0032 via a bisection 
search (and with a fractional precision of ?s 10~ 15 ). The other 
parameters defining the gaussian initial profile of the massless 
field are ro = 40 and a = 5. The snapshots show that the real 
scalar field enters the region containing the bulk of boson star 
at t ~ 22, implodes through the origin at t m 45, leaves the 
boson star region at t ~ 70, and, finally, completely disperses 
from the computational domain at t ~ 100. The boson star 
enters the critical state at roughly the same time that the real 
field leaves the domain, and remains in that state for a period 
of time which is long compared to the crossing time of the 
massless field. At t ta 510, the boson star begins to depart 
significantly from the critical state. 



hole formation. 



C. Final Fate of Subcritical Evolutions 

In previous work on the problem of critically perturbed 
spherically symmetric boson stars [HQ, it was conjec- 



FIG. 2: Time evolution of the central value of the scalar field 
modulus for subcritical evolution of perturbed boson stars. 
The figure shows the time evolution of \4>(t, 0)| for marginally 
subcritical evolutions generated from boson star initial states 
with 0o (0) = 0.035, 0.04 and 0.05. See the text for a descrip- 
tion of key features of this plot. 



tured that the end state of subcritical evolution was char- 
acterized by dispersal of the boson star to large distances 
(relative to the size of the initial, stable star). This con- 
jecture was at least partially influenced by the behaviour 
observed, for example, in the collapse of a massless scalar 
field 3 , where subcritical evolutions do involve complete 
dispersal of the field. However, another key reason for 
what we claim is a misidentification of the true subcriti- 
cal end-state, was that the simulations described in [l|, Q 
simply were not carried out for sufficient coordinate time 
to observe the nature of the late-time dynamics. Our 
current simulations strongly suggest that subcritical evo- 
lutions lead to a "relaxation" of the critically perturbed 
state to something that approximates a boson star (not 
necessarily the original star) undergoing large amplitude 
oscillations. As argued in the next subsection, these os- 
cillations can largely be identified with the fundamental 
perturbative mode associated with the final boson star 
state. The numerical evidence also suggests that, at least 
in many cases, these oscillating configurations eventually 
re-collapse and form black holes; a "prompt" re-collapse 
can be seen in the <?!>o(0) = 0.05 data in Fig. 

Fig. [5] displays the long-time behaviour of 
max r (2M(t, r)/r), \<f>(t, 0)| and ip(t,0) for a near- 
critically perturbed boson star (0o(O) = 0.04,^3 « 
0.00342) for r max = 200 (with mesh spacing 
Ar = 200/4096 w 0.049). Note that this is a sub- 
critical evolution, so that a black hole does not form. 
As shown in more detail in previous figures, the boson 
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FIG. 3: Transition of perturbed boson stars in critical evolu- 
tions. The solid curve shows the parametric mass vs radius 
plot of static boson stars (curve parameter 0o(O) increasing 
from right to left), where we have defined the stellar radius, 
R, so that M(R) = 0.99 M(oo) = 0.99 Madm- Triangles la- 
bel the initial configurations, squares show the corresponding 
critical solutions (identified as one-mode-unstable boson stars 
with oscillations that are largely in the fundamental mode), 
and the dashed arrows represent schematically the transition 
between the initial and critical states. See the text for more 
details. 



FIG. 4: Measured lifetime scaling laws for critically perturbed 
boson stars. This figure shows the measured lifetimes of vari- 
ous near-critical evolutions of perturbed boson stars as a func- 
tion of In |A 3 - A|], for cases with <j)(0) = 0.02, 0.035, 0.04 and 

0. 05. Quoted scaling exponents, 7 (see (|67[) ). are computed 
from linear least-squares fits to the data. The apparent con- 
vergence of the data for different (j>o(0) as In \A$ — AJ| — > is 
not significant, as it reflects calculations far from criticality 

1. e. far from the In \ A$ — Af| — > — 00 limit. See the text for 
additional details. 



star enters a critical state (well approximated by an 
unstable boson star) shortly after the real scalar field 
leaves the computational domain (t « 100). While in the 
critical state, the star oscillates with what we assume is 
the frequency of the first harmonic, as computed from 
perturbation theory using, the unstable boson star state 
as the background (see [3, [1]). At t « 300 the star 
begins to evolve away from the more compact critical 
configuration, decreases in central density, expands in 
size, and starts to pulsate with a different frequency. 
Although at late time the oscillation amplitudes are 
much larger than those seen in the critical phase of 
evolution, we will show in the following section that 
the oscillations can nonetheless be largely attributed 
to excitations of the fundamental perturbative mode 
associated with the final boson star state. 

Fig. [5] shows the long-time behaviour of the modulus 
of the central value of scalar field, |</>(i,0)|, for initial 
configurations with (j>o(0) = 0.035,0.04 and 0.05, with 
''max = 100, but with Ar maintained at 50/1024 as in 
Fig. [5l Again, we use A 3 to tune the evolution of the bo- 
son stars to criticality and the figure shows a marginally 
subcritical evolution. In general, the computed value of 
A3 is a function of r max , as is the specific stable boson 
star to which the critical evolution relaxes. However, the 



results shown in the figure support our claim that an 
oscillatory phase (rather than dispersal) generically fol- 
lows near-critical evolution of driven boson stars in the 
marginally subcritical case. 

Fig.[7]shows the long time behaviour of subcritical evo- 
lution of the modulus of the central scalar field value, 
\4>(t, 0)|, with an initial boson star given by </>o(0) = 0.04. 
Here, we vary the position of the outer edge of the 
computational domain r max , while keeping the resolu- 
tion, Ar, fixed at 50/1024 as previously. For each of 
r ma i = 50, 100, 200 and 400, we tune A 3 to generate 
a critical evolution (the specific values of A3 obtained 
are listed in Table I). This set of calculations provides 
evidence for the convergence of the critical solution (in- 
cluding the critical value of the control parameter, A3), 
as r max — > 00 at fixed resolution. This in turn strongly 
suggests that the final oscillatory states identified in sub- 
critical evolutions are not artifacts of our use of a finite 
computational domain. 

In order to illuminate the nature of typical post-critical 
oscillations, Fig. [5] shows the square of the discrete fast 
Fourier transform, T [\<fi(t, 0)|], of the central scalar field 
modulus for the same set of simulations used to prepare 
Fig. [7] The transform is taken for discrete times, t n , 
satisfying 2500 < i" < 7000, a period when the boson 
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FIG. 5: Long time behaviour of subcritical evolution for 
0(0, 0) = 0.04 with r max = 200. This figure shows the long- 
time behaviour of max r (2M(t, r)/r), |(/>(t,Q)| and ij)(t,0) for 
a near-critically perturbed boson star (</>o(0) = 0.04,^3 » 
0.00342). The left side of the figure shows the evolution of 
the perturbed star in its critical state (100 < t < 300), and 
the evolution shortly after the star leaves its critical state. 
The right side of the figure focuses on oscillations seen at 
later times 1000 < t < 7680. These plots provides evi- 
dence that the final state of subcritical evolution is charac- 
terized by large amplitude oscillations about something ap- 
proximating a boson star on the stable branch, rather than 
dispersal of the complex field as suggested in 1, 2]. De- 
tailed calculation (see Sec. MI D[) shows that the pulsation 
frequency is approximately the fundamental mode frequency 
computed from perturbation theory about a background sta- 
ble boson star solution with 4>o(0) = 0.023. We also note 
the overall lower-frequency modulation of the post-critical os- 
cillations. This effect is not yet understood, although one 
possible explanation — namely that the envelope modulation 
represents "beating" of the fundamental and first harmonic 
modes — appears to be ruled out. 



star has undergone the transition from critical evolution 
to post-critical oscillation. The figure clearly shows the 
convergence of the fundamental mode oscillation, as well 
as a first harmonic. The next section provides a more 
detailed analysis of the observed fundamental mode ex- 
citations. 



D. Perturbation Analysis of Subcritical Oscillations 

We now proceed to an application of perturbation 
theory to the oscillations seen in long-time evolutions 
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FIG. 6: Long time behaviour of subcritical evolution with ini- 
tial configurations <^o(0) = 0.035,0.04 and 0.05, for r max = 
100. The figures show the modulus of the central scalar 
field values, [</>(£, 0)|, vs time, using the same resolution 
Ar = 50/1024 used to generate the data shown in Fig. [5] 
Each of the three distinct boson stars is driven to a different 
critical solution, and subsequently relaxes to a different final 
oscillatory state. This provides evidence that the final end 
state of marginally subcritical evolution in generic driven bo- 
son stars does not involve dispersal of the bulk of the complex 
field to infinity. 



of marginally subcritical configurations, such as those 
shown in Fig. Here we follow HI and 0, and refer 
the interested readers to those sources for details of the 
approach that we do no include here. In particular, we 
emphasize that we have not carried out the complete per- 
turbation analysis ourselves, but are simply using a com- 
puter code provided by Hawley [f| to analyze our current 
simulations. Nonetheless, to make contact between the 
perturbative and simulation results, it is useful to briefly 
review the setup of the perturbative problem. 

To formulate the equations for the perturbation anal- 
ysis, we first rewrite the complex scalar field as 

<P(t,r) = (Mt,r)+iMt,r))e-^\ (68) 

(Note that this representation is distinct from <p = 
4>i + i4>2, and the reader should be careful not to con- 
fuse the ip' s used here with the conformal metric vari- 
able, ip.) Additionally, the spacetime metric is written in 
Schwarzschild-like (polar-areal) coordinates: 

ds 2 = -e^dt 2 + e X{t ^dr 2 

+ r 2 (d6 2 + sin 2 6d(p 2 ) . (69) 

We further introduce four perturbation fields, 
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FIG. 7: Long time behaviour of subcritical evolution with an 
initial boson star characterized by <po(0) = 0.04, for r max = 
50, 100, 200 and 400. The figures show the modulus of the cen- 
tral scalar field values, |$>(t, 0)|, vs time, with the resolution 
Ar fixed at 50/1024 as in previous figures. The evolutions 
are tuned to criticality for different r max (see Table I). The 
top figure shows the overall evolutions for r max = 50, 100, 200 
and 400 from t = to t = 7680. The middle figure focuses on 
the evolution of the perturbed boson star during the period of 
near-critical evolution, 70 < t < 270, for the cases r max = 200 
and 400. The near coincidence of the two curves in this case 
provides strong evidence for convergence of our calculations 
(at fixed spatial resolution) as r max —+ oo . The bottom figure 
focuses on the late time evolution — 200 < t < 7680 — again for 
r max = 200 and 400, and provides additional support for our 
claim that the final oscillatory states we observe in subcritical 
evolution are not an artifact of the use of a finite computa- 
tional domain. 



FIG. 8: Long time behaviour of subcritical evolution with 
an initial boson star characterized by 4>o(0) = 0.04, for 
rmax = 50, 100, 200 and 400. The figures show the square 
of the (discrete) Fourier transform T [\(j>(t, 0)|], of the central 
scalar field modulus, using the same calculations described 
in Fig. [7] The transform is taken from a data set defined 
at 691 discrete times, t n satisfying 2500 < t n < 7700, dur- 
ing which time the critically perturbed boson star is in its 
final oscillatory state. Again, the resolution, Ar = 50/1024, 
is the same used in previous calculations The fundamental 
mode computed for the case r max = 200 is approximately 
uj ~ 33 x 6 x 10~ 4 = 0.0198, in good agreement with our 
perturbation-theory estimate computed in Sec. IIII Dl The 
figure shows that the computed frequency of the fundamental 
mode converges for increasing r max . The graph also shows 
evidence for at least one higher overtone which persists as 
r max — » oo. The figure inset shows the overall amplitudes of 
the computed Fourier components. 



SX(t, r), 5v{t, r), Si^i (i, r) and Sip-i{t,r), which repre- 
sent the perturbations about the equilibrium values 
Xo{r),v a {r),(j) {r): 

A(t,r) = X (r)+5X(t,r), (70) 

u(t,r) = v {r) + Sv(t,r) , (71) 

Mt,r) = Mr)(l + $Mt,r)) , (72) 

V> 2 (t,r) = Mr)5Mt,r). (73) 



With the above definitions we can write the coupled 
Einstein-Klcin-Gordon field equations as a set of PDEs 
for the functions SX, 5v, Sipi and <5"02- With some manip- 
ulation we can then eliminate 8v and Sip2 to produce a 
system of two coupled second-order PDEs for 6ipi and 
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Note that these equations involve only second time 
derivatives (i.e. there are no terms involving Stpi or SX), 
and that they are linear in the second time derivatives. 
If we thus assume a harmonic time-dependence for the 
perturbed fields: 



5Xi{t,r) 



<5Vi(r)e i(Tt 
5Xi{r)e iat . 



(76) 
(77) 



then the equations for the perturbations contain a only 
in the form a 2 , and the sign of a 2 , as computed by solv- 
ing a particular mode equation, determines the stability 
of that mode. (Note that the system can be shown to 
be self- adjoint so that the values of a 2 must be real.) 
If any of the values of a 2 are found to be negative, then 
the associated perturbations will grow and the boson star 
will be unstable. Moreover, as the eigenvalues form an 
infinite discrete ordered sequence, examining the funda- 
mental radial mode <jq 2 determines the overall stability of 
any particular star with respect to radial perturbations. 

In order to compare the simulation results with those 
given by perturbation theory, we first observe that there 
is a difference in the choice of the time coordinates used 
in the two calculations. Specifically, in the perturbative 
analysis 0, [H|, the lapse is chosen to be unity at the 
origin, so we have 



perturbative 



simulation 



We also note that there is a factor of 2 difference in the 
definitions of T M „ used in the two calculations, and that 
the definition of the complex field, (f)(t,r), in the per- 
turbative calculation includes a factor of v8tt- We thus 
have 



perturbative 



47T0 



simulation 



The numerical technique for obtaining the fundamen- 
tal mode and first harmonic mode frequencies of boson 
stars has already been described in Q and will not be 
repeated here; again, we will simply quote and use re- 
sults from that study. From Fig. [5] we note that there 



are 10 oscillations between t = 2553.8 and t = 5583.8, 
giving a period T ~ 333. Hence we have an oscilla- 
tion frequency a = 2tt/T sw 0.019. The time average 
of the lapse function, (a(t,Q))t, in the interval is 0.89, 
and so a 2 jo? w 0.00045. We also compute the time 
average of <p(t, 0) in the interval, and use the resulting 
value to identify the stable boson star solution about 
which we perform the perturbation analysis. We find 
(0o(*,O)) t w 0.023 x Vfe = 0.0815. For a boson star 
with 4>o(0) = 0.0815, the perturbative calculations (see 
Fig. 7 of Q) predict a 2 = 0.00047, which is in reasonable 
agreement with the simulation results. Hence the oscil- 
lations that occur in the post-critical regime appear to 
be largely fundamental mode oscillations of a final-state, 
stable, boson star. We also remark that since the oscilla- 
tions are of such large amplitude, it does not appear pos- 
sible to precisely identify an effective background state 
(i.e. an effective value of <^>o(0)), so the level of agreement 
in the oscillation frequencies is probably as good as one 
could expect. 



IV. SUMMARY 

We have investigated type I critical phenomena of 
ground state boson stars in maximal-isotropic coordi- 
nates by perturbing the stars with in-going pulses of 
a real scalar field. In particular, contrary to previous 
claims, we find that the end state of generic subcritical 
evolution is a stable boson star executing large amplitude 
oscillations, and that the oscillations can be largely un- 
derstood as excitations of the fundamental normal mode 
of the end-state star. For the particular example that we 
examined in detail, the oscillation frequency of the post- 
critical state was estimated to be a 2 /a 2 w 0.00045, in 
good agreement with the frequency of the fundamental 
mode computed in perturbation theory, <7g = 0.00047. 
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APPENDIX A: FINITE DIFFERENCE 
ALGORITHM 

Here we present the details of the numerical method 
used in our computations. We solve the PDEs (fT2"|) - 
(p~9j) by a finite difference method. We replace the 
(t, r) continuum by a discrete lattice of grid points, 
and approximate the continuum field quantities J- = 
{a,P,^p,K r r ,(j)i,^i,Ui}, where i = 1,2,3, by a set of 
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FIG. 9: Stencil for an 0(h 2 ) Crank-Nicholson scheme for a 
PDE in one space dimension and time. 



grid functions T h = {a h , (3>\ ip h ,K r ^, $, $f , U^} which 
are solutions of the finite difference approximation (FDA) 
of the PDEs. Denoting the uniform (constant) spa- 
tial and temporal mesh spacings by Ar and At, respec- 
tively, the finite difference grid is given by (t n , Tj ) where 
r.j = r + {j - l)Ar, j = 1, ••• ,N r and t n = nAt, 
n = 0, • • ■ ,Nt- For any grid function u h £ T^, the value 
at (t n ,rj) is denoted by it™ and is an approximation of 
the continuum value u(t n ,rj). 

In discretizing evolution equations ([14")) — (fTo]) we make 
exclusive use of Crank-Nicholson schemes, with second 
order spatial differences. The key idea of a Crank- 
Nicholson method is to keep the differencing centred in 
time as well as in space, and a typical stencil used for such 
a scheme is illustrated in Fig. [9] The constraint equations 
(|12p and (|13p are coupled, nonlinear, ordinary differen- 
tial equations, and following 0(Ar 2 ) finite differencing 
(see below) and are solved using a point-wise Newton's 
method. That is, at each grid point, (t n ,rj), we solve 
for the pair (-0™, (K r r )™) using Newton's method for two 
equations in two unknowns. The slicing condition (jTSj) 
is linear, so, after being discretized using second-order 
finite differences, can be solved directly using a tridiag- 
onal solver. Finally, once the values a™ and K r r ™ have 
been computed, an 0(Ar 2 ) discretization of §T§\i is easily 
integrated to yield the /?". 

To aid in the presentation of the finite difference equa- 
tions, it is convenient to define the following difference 
operators: 
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and the averaging operator 
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We also define /j ± , which has the same definition as /i ± , 
but which has a higher precedence over other algebraic 
operations, e.g., 
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The FDAs of the Klein-Gordon equations can then be 
written as: 



+ j" 



n+l 



At 



A + (&). = M+ (-H, • ^ 



3 



(Al) 



A n w. 

o j 



j+i j-i 
2Ar 



A + ^);;=^ A oi^+^ n 



(A2) 



15 



. t , .71 t 

a, (n,), = P-i 



A„ 



1 (v> 4 )™ 



rV ( /?n, + — $, 



Q' 



[a^ 2 mVi(l-&)]' 



(aJT r ) . + 2/3 



3 ' 3 (rib 2 ) n 



(Hi) . (A3) 



4(0i)J-(0i)a\ 


= o, 


(All) 


. n 


= o, 


(A12) 


4(n,-)J - (n,-)S 


= o, 


(A13) 



/ \ 



for all i and n. The outer boundary conditions are 



where i = 1, 2, 3. 

The FDA of the Hamiltonian constraint is 



<J> 

£ n £ / . r n * . 

4 3 + \ 06 J r . 



= 0, (A14) 



3 



V j o *j / 16 v y 3 

/E- =1 (^ + n?) , , 



V> 4 



i=l 



and the FDA of the momentum constraint is 

fcjii. f A_(K r r ). + 3A_(r^ 2 ) j fl_ i^—^j 



(A5) 
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Similarly, the FDAs for the maximal-isotropic conditions 
are 
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We also adopt a scheme for numerical dissipation given 
by Kreiss and Oliger In other words an additional 
term 
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' + V +KO l j 



is added to the right hand side of (IA2|) , for 3 < j < N r — 2 
(and similarly to the right hand side of (|A3|) for the IT), 

where A +KQ is defined by 
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Here, is an adjustable parameter satisfying < td < 1, 
and is typically chosen to be 0.5. We note that the addi- 
tion of Kreiss-Oliger dissipation changes the truncation 
error of the FDAs at 0(At 3 , Ar 3 ) and thus does not effect 
the leading order error of a second order (0(Ai 2 , Ar 2 )) 
scheme. The dissipation is useful for damping high fre- 
quency solution components that are often associated 
with numerical instability. 
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APPENDIX B: CONVERGENCE TESTING 
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respectively, where r. l = (rj + fj-_i)/2. 
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The regularity conditions are implemented as 
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Here we present the results of a convergence test of the 
code that evolves boson stars in spherical symmetry. 

In Fig.[Tn]we plot the mass aspect function at the outer 
boundary of the computational domain, M(i,r max ), as a 
function of time, and from four computations with grid 
spacings, Ar, in a 8:4:2:1 ratio. As was the case for the 
calculations discussed in the main text Sec. IIII1 our con- 
vergence study uses a pulse of massless scalar field im- 
ploding onto a stable boson star So long as no scalar field 
(either real or complex) propagates off the computational 
grid, M(t, r max ) should be constant in time (and equal to 
the ADM mass), in the limit that Ar -> (with At — > 
implied since A is always held fixed as Ar is varied). 
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In our test, the boson star has a central field value, 
4>o = 0.01, while the incoming massless scalar field pulse 
is a gaussian of the form (j6"5j) with A 3 — 0.001, r = 40 
and a = 3. The outer boundary is r max = 300, and 
we compute with N r = 1025,2049,4097 and 8193. Dur- 
ing the time interval 40 < t < 50, the real scalar 
field is concentrated near the origin and interacts most 
strongly with the complex field. This results in a lo- 
calized fluctuation of the computed ADM mass that is 
evident in the plots. However, M(t, r max ) clearly tends 
to a constant value as the resolution is increased. In ad- 
dition, from the differences of -M(i,r max ) computed at 
different resolutions (e.g. M Ar (t, r max ) - M 2Ar {t, r max ), 
M 2Ar (t, r max ) - M 4Ar (t, r max ), etc.), we find strong evi- 
dence that the overall difference scheme is converging in 
a second order fashion. 
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FIG. 10: Convergence test of the spherically symmetric code. 
The estimated ADM mass, M(t,r max ), is plotted against 
time, t, for four calculations using numbers of spatial grid 
points, N r , of 1025,2049,4097 and 8193, so that the corre- 
sponding mesh spacings, Ar, are in a 8:4:2:1 ratio. The ini- 
tial data parameters for the computations are: (j>o = 0.01 for 
the complex field, and A3 — 0.001, ro = 40 and a — 3 for 
the massless field (see (|65jl ) . The mass decreases with time in 
general, with a significant fluctuation at 40 < t < 50, when 
the real scalar field is close to the origin and strongly inter- 
acts with the boson star. The variation in the computed total 
mass tends to vanish as we go to higher resolution. Combin- 
ing results from the four calculations we find strong evidence 
that the finite difference scheme is second order accurate as 
expected. 
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